Diverse Phases of Carbonaceous Materials from Stochastic Simulations

Amorphous carbon systems are emerging to have unparalleled properties at multiple length scales, making them the preferred choice for creating advanced materials in many sectors, but the lack of long-range order makes it difficult to establish structure/property relationships. We propose an original computational approach to predict the morphology of carbonaceous materials for arbitrary densities that we apply here to graphitic phases at low densities from 1.15 to 0.16 g/cm3, including glassy carbon. This approach, dynamic reactive massaging of the potential energy surface (DynReaxMas), uses the ReaxFF reactive force field in a simulation protocol that combines potential energy surface (PES) transformations with global optimization within a multidescriptor representation. DynReaxMas enables the simulation of materials synthesis at temperatures close to experiment to correctly capture the interplay of activated vs entropic processes and the resulting phase morphology. We then show that DynReaxMas efficiently and semiautomatically produces atomistic configurations that span wide relevant regions of the PES at modest computational costs. Indeed, we find a variety of distinct phases at the same density, and we illustrate the evolution of competing phases as a function of density ranging from uniform vs bimodal distributions of pore sizes at higher and intermediate density (1.15 g/cm3 and 0.50 g/cm3) to agglomerated vs sparse morphologies, further partitioned into boxed vs hollow fibrillar morphologies, at lower density (0.16 g/cm3). Our observations of diverse phases at the same density agree with experiment. Some of our identified phases provide descriptors consistent with available experimental data on local density, pore sizes, and HRTEM images, showing that DynReaxMas provides a systematic classification of the complex field of amorphous carbonaceous materials that can provide 3D structures to interpret experimental observations.

The root mean square deviation (RMSD) measures the deviation of a structure from a reference structure (RMSD=0.0 indicates a perfect overlap). The RMSD is defined as: where N is the number of atoms, m i is the mass of atom i, X i is the coordinate vector for target atom i, Y i is the coordinate vector for reference atom i, and M is the total mass. If the RMSD is not massweighted, all m i = 1 and M = N.
The root mean square fluctuation (RMSF) of a structure is the time average of the deviation of a structure from its average value. It is calculated according to the equation below, where x i is the coordinates of particle i, and ⟨x i ⟩ is the ensemble average position of i: The RMSD quantifies how much a structure diverges from a reference over time, the RSMF reveals which areas (atoms) of the system are the most mobile. While RMSD is frequently calculated with respect to the initial state, the RMSF is calculated with respect to an average structure of the simulation. An area of the structure with high RMSF values frequently diverges from the average, indicating high mobility. Note that we compute RMSD and RMSF after a rigid-body optimal alignment of the atomic coordinates in each snapshot to a reference structure (starting or average configuration).

Definition of massaging protocol
The parameters we have chosen to explore in the 'massaging' procedure, reported in Table 1, belong to the atom, bond, valence angle and torsion sections of the ReaxFF force field, which is divided into 7 parts, containing the general, atom, bond, off-diagonal, valence angle, torsion angle, and hydrogen bond parameters -we refer to the ReaxFF manual 3,4,5 for a detailed description. As explained in the ReaxFF manual, all parameters without physical meaning are named after the partial energy contribution (valp1, vapl2, etc.). In contrast, the other ones, like the torsional rotational barriers (V1, V2, V3), have names connected to physical quantities or descriptors.

Figure S2
. Initial and final structures obtained after 20 ps massaging MD at T=2000 K in the NVT ensemble. The massaged parameters are indicated with P i in each case, and the molecular structures are represented with cyan lines connecting the carbon atoms and filled colored rings. Different colors are used to distinguish the number of ring members. The initial structure is the same prototype system, containing 4176 atoms (simulation box = 42 × 38 × 45 Å 3 ), cited in Figure S1. Table S2. Force field parameters randomly selected for the four-step-massage NVT MD simulations at T=2000 K (50% parameter reduction), starting from the small prototype model of a carbonaceous material. Figure S3. Root mean square deviation of all the atoms of the configurations sampled at the end of each massage (forty in all) relative to the initial structure. The initial structure is the same prototype system containing 4176 atoms (simulation box = 42 × 38 × 45 Å 3 ) shown in Figure S2. Each NVT MD simulation was carried out at 2000 K for 50 ps. The modified parameters are reported in Table  S2. Table S3. Structural descriptors of the geometries of the prototype structure obtained at the end of the 4-step-massage and after the GO procedure. D=Box density (g/cm 3 ) C 1 =Fraction of C atoms with coordination 1. C 2 =Fraction of C atoms with coordination 2. C 3 =Fraction of C atoms with coordination 3. C 4 =Fraction of C atoms with coordination 4. A (sp) =Average difference between ideal 180° angle for sp C atoms and effective angles formed by those C atoms characterized by coordination 2. Large values correspond to "bended wires" deviating from ideal linearity (degrees). A (sp2) =Average difference between the ideal 120° angle for sp 2 C atoms and effective angles formed by those C atoms characterized by coordination 3. Large values correspond to "bended sheets" deviating from perfect planarity of ideal graphitic foils (degrees). R 5 =Ratio between the number 5-member rings and the total number of atoms R 6 = Ratio between the number 6-member rings and the total number of atoms R 7 = Ratio between the number 7-member rings and the total number of atoms CN 5 =coordination number (second sphere) sasa=Solvent-accessible surface area (Å 2 ) D ave =Average Density of the System (g/cm 3 ) FV=Fraction of Free-Volume (%) Figure S4. Final representative structures of the ten four-step-massages described in Tables S2-S3. The initial structure is the same prototype system containing 4176 atoms (simulation box = 42 × 38 × 45 Å 3 ) shown in Figure S2. PLD = Pore limiting diameter, Å LCD = Largest cavity diameter, Å S AC,T = Total accessible surface area, m 2 /g S AC,A = Network accessible surface area, m 2 /g V PO,T = Total probe-occupiable volume, cm 3 /g V PO,A = Network accessible probe-occupiable volume, cm 3 /g F He,T = Helium pore volume fraction C 1 =Fraction of C atoms with coordination 1. C 2 =Fraction of C atoms with coordination 2. C 3 =Fraction of C atoms with coordination 3. C 4 =Fraction of C atoms with coordination 4. A (sp) =Average difference between ideal 180° angle for sp C atoms and effective angles formed by those C atoms characterized by coordination 2. Large values correspond to "bended wires" deviating from ideal linearity (degrees). A (sp2) =Average difference between the ideal 120° angle for sp 2 C atoms and effective angles formed by those C atoms characterized by coordination 3. Large values correspond to "bended sheets" deviating from perfect planarity of ideal graphitic foils (degrees). R 5 = 100 x Ratio between the number 5-member rings and the total number of atoms R 6 = 100 x Ratio between the number 6-member rings and the total number of atoms R 7 = 100 x Ratio between the number 7-member rings and the total number of atoms CN 5 =coordination number (second sphere) sasa=Solvent-accessible surface area (Å 2 ) D ave =Average Density of the System (g/cm 3 ) FV=Fraction of Free-Volume (%)  3 mass density (MM1/MM4 is simplified to m1m4 in the inset for the sake of notation, and analogously for the other combinations). PLD = Pore limiting diameter, Å LCD = Largest cavity diameter, Å S AC,T = Total accessible surface area, m 2 /g S AC,A = Network accessible surface area, m 2 /g V PO,T = Total probe-occupiable volume, cm 3 /g V PO,A = Network accessible probe-occupiable volume, cm 3 /g F He,T = Helium pore volume fraction C 1 =Fraction of C atoms with coordination 1. C 2 =Fraction of C atoms with coordination 2. C 3 =Fraction of C atoms with coordination 3. C 4 =Fraction of C atoms with coordination 4. A (sp) =Average difference between ideal 180° angle for sp C atoms and effective angles formed by those C atoms characterized by coordination 2. Large values correspond to "bended wires" deviating from ideal linearity (degrees). A (sp2) =Average difference between the ideal 120° angle for sp 2 C atoms and effective angles formed by those C atoms characterized by coordination 3. Large values correspond to "bended sheets" deviating from perfect planarity of ideal graphitic foils (degrees). R 5 = 100 x Ratio between the number 5-member rings and the total number of atoms R 6 = 100 x Ratio between the number 6-member rings and the total number of atoms R 7 = 100 x Ratio between the number 7-member rings and the total number of atoms CN 5 =coordination number (second sphere) sasa=Solvent-accessible surface area (Å 2 ) D ave =Average Density of the System (g/cm 3 ) FV=Fraction of Free-Volume (%)  3 mass density (MM1/MM4 is simplified to m1m4 in the inset for the sake of notation, and analogously for the other combinations). PLD = Pore limiting diameter, Å LCD = Largest cavity diameter, Å S AC,T = Total accessible surface area, m 2 /g S AC,A = Network accessible surface area, m 2 /g V PO,T = Total probe-occupiable volume, cm 3 /g V PO,A = Network accessible probe-occupiable volume, cm 3 /g F He,T = Helium pore volume fraction C 1 =Fraction of C atoms with coordination 1. C 2 =Fraction of C atoms with coordination 2. C 3 =Fraction of C atoms with coordination 3. C 4 =Fraction of C atoms with coordination 4. A (sp) =Average difference between ideal 180° angle for sp C atoms and effective angles formed by those C atoms characterized by coordination 2. Large values correspond to "bended wires" deviating from ideal linearity (degrees). A (sp2) =Average difference between the ideal 120° angle for sp 2 C atoms and effective angles formed by those C atoms characterized by coordination 3. Large values correspond to "bended sheets" deviating from perfect planarity of ideal graphitic foils (degrees). R 5 = 100 x Ratio between the number 5-member rings and the total number of atoms R 6 = 100 x Ratio between the number 6-member rings and the total number of atoms R 7 = 100 x Ratio between the number 7-member rings and the total number of atoms CN 5 =coordination number (second sphere) sasa=Solvent-accessible surface area (Å 2 ) D ave =Average Density of the System (g/cm 3 ) FV=Fraction of Free-Volume (%)

Convergence Check
For the three examined densities (1.15, 0.50, and 0.16 g/cm 3 ) and for all combinations of massages {MM1,MM2,MM3} x {MM4,MM6,MM8}, the final structure of the last massage was equilibrated for about fifty picoseconds (ps) using the original (unbiased, not-massaged) ReaxFF (C.ff) force field at T=2000 K. The average configuration obtained from the final ten ps of this equilibration was used for the subsequent GO steps whence the final descriptor analysis. To check convergence, the root mean squared deviations (RMSD) of all the atoms from the average structure were evaluated in the last ten ps of the equilibration dynamics: these are illustrated in Figure S5(left panel) focusing on the intermediate density (0.50 g/cm 3 ). A superimposition of the average structure with another structure having a RMSD of approximately 7.0 Å obtained during the last ten ps of the dynamics is also shown in Figure S5(right panel). Inspection of Figure S5 indicates that the RMSD are small and oscillate in a range of 6.2-8.5 Å. Visual inspection of the superimposition further confirms the validity of convergence. Figure S5. On the left, distribution of the RMSD of all the structures sampled during the last ten ps of the equilibration dynamics with the C.ff force field at T=2000K relative to the average configuration. Results refer to mass density of 0.50 g/cm 3 and are reported for all combinations of massages: {MM1,MM2,MM3} x {MM4,MM6,MM8} (MM1/MM4 is simplified to m1m4 in the inset for the sake of notation, and analogously for the other combinations). On the right, a typical superimposition between all the atoms of the average structure and those of a configuration sampled during the last ten picoseconds of the equilibration dynamics having a RMSD=7 Å.                Figure S22. A narrow tunnel identified by the CAVER software in one of the sampled structures exhibiting a bimodal distribution of pore sizes. The tunnel corresponds to the continuous green line, and for the sake of clarity, we have highlighted with a pink aura its intersection with the wall of the periodic unit cell.